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1.  INTRODUCTION 

Accurate  modeling  of  tropical  cyclone  motion  and  in¬ 
tensity  change  requires  resolving  the  flow  both  within  and 
around  the  storm.  Since  the  spatial  scales  in  these  two 
regions  differ  substantially,  uniform  resolution  is  inherently 
inefficient:  the  grid  should  be  refined  only  near  the  storm. 
This  idea  motivates  conventional  nested-grid  methods  such 
as  used  in  VICBAR  (DeMaria  et  aL,  1992)  and  the  GFDL 
model  (Kurihara  et  aL,  1998),  Adaptive  multigrid  methods 
also  achieve  nonuniform  resolution  by  superimposing  uni¬ 
form  grids  of  different  mesh  sizes,  but  they  combine  this 
idea  with  multigrid  processing  (Brandt,  1977)  to  achieve 
optimum  solution  speed  and  provide  accurate  truncation 
error  estimates.  The  latter  can  be  used  in  an  adaptive  mesh 
refinement  scheme  to  provide  just  the  resolution  needed  at 
each  point. 

The  MUDBAR  model  of  Fulton  (2000)  demonstrates 
the  potential  of  adaptive  multigrid  methods  in  the  context 
of  a  nondivergent  barotropic  model.  We  now  consider  the 
extension  of  these  techniques  to  the  next  level  of  dynamical 
complexity,  i.e.,  the  shallow-water  equations. 

2.  MODEL  DESCRIPTION 

In  this  section  we  briefly  describe  the  model;  further 
details  can  be  found  in  Mitchell  and  Fulton  (2000). 

2.1  Governing  Equations 

We  work  with  the  shallow-water  equations  in  the  form 

Ut  +  UUx  +  VUy  -  fv  +  (px  =  0, 

Vt  +  UVx  +  VVy  -I-  /w  -{-  =  0,  (1) 

+  (<J>  +  (t>){Ux  +  Vy)  +  U<l>x  +  V(j>y  =  0, 

where  u  and  v  are  the  velocity  components  in  x  and  y, 
respectively,  <f>  is  the  deviation  of  the  geopotential  gh  from 
a  constant  reference  value  <t>,  and  /  is  the  Coriolis  param¬ 
eter.  The  model  domain  is  a  rectangle  on  a  /3-plane.  At 
the  boundaries  we  use  open  boundary  conditions,  specify¬ 
ing  the  quantity  Vn  —<p/c  (where  Vn  is  the  outward  normal 
velocity  component  and  c  =  y/^)\  the  tangential  velocity 
component  is  also  specified  where  there  is  Inflow. 


at  a  new  time  level  take  the  form 

u  -h  M<j)x  =  F, 

+  M<t>y  =  G,  (2) 

(p  +  <t>At{Ux  Vy)  =  H, 


where  At  is  the  time  step  and  the  right-hand  side  depends 
on  values  at  the  previous  time  levels.  The  initial  time  step 
(second-order  Runge-Kutta)  takes  a  similar  form.  Eliminat¬ 
ing  u  and  V  leads  to  a  Helmholtz  problem 

<t>  -  (cAt)^(0xx  +  <l>yy)  =  Q  =  H-  +  Gy)  (3) 

to  be  solved  for  the  geopotential;  the  corresponding  velocity 
components  can  then  be  obtained  from  (2). 

2.3  Space  Discretiza tion 

The  model  is  discretized  in  space  using  second-order 
centered  finite  differences  on  an  Arakawa  C-grid  with  mesh 
spacing  h.  The  discretization  of  (2)  is  straightforward  (but 
requires  one-sided  differences  on  the  right-hand  side  when 
discretized  at  boundary  points).  The  discretization  of  (3) 
takes  the  form 

=  (4) 

where  the  discrete  operator  has  the  difference  stencil 
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(5) 


with  7  =  cAt/h  being  the  Courant  number. 

To  apply  the  boundary  conditions  we  average  the  nor¬ 
mal  velocity  across  the  boundary  (using  a  ghost  point  value) 
and  eliminate  the  velocity  components  using  the  discretized 
momentum  equations  to  obtain  a  modified  form  of  (4)  ap¬ 
plied  at  the  boundary.  For  example,  (5)  is  replaced  by 
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(6) 


2.2  Time  Discretization 

To  permit  using  larger  time  steps  than  would  be  al¬ 
lowed  by  the  CFL  condition  of  an  explicit  method,  we  use 
a  semi-implicit  (leapfrog/trapezoidal)  time  discretization. 
The  resulting  implicit  equations  for  the  solution  (u,v,0) 
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at  the  northeast  corner;  the  right-hand  side  is  modified 
similarly.  Note  that  while  the  system  (2)  has  a  staggered 
discretization,  the  discrete  Helmholtz  problem  (4)  involves 
only  values  at  the  gridpoints. 
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2.4  Multigrid  Solution 

To  solve  (4)  efficiently  we  use  a  multigrid  method. 
This  uses  Gauss-Seidel  relaxation,  full  weighting  of  residu¬ 
als,  and  bilinear  interpolation  of  corrections  in  a  Full  Multi¬ 
grid  (FMG)  algorithm  with  bicubic  initial  Interpolation. 
Since  the  boundary  conditions  are  combined  with  the  inte¬ 
rior  equations  as  described  above,  proper  residual  transfers 
must  be  used  at  the  boundaries.  For  example,  the  standard 
full  weighting  stencil 
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is  replaced  by 
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at  the  west  boundary  and  northwest  corner,  respectively. 
For  large  time  steps  (7  ^  1)  the  operator  approaches 
the  Laplacian  operator,  and  the  multigrid  smoothing  factor 
approaches  Jx  =  0.25  (with  red-black  ordering);  for  smaller 
time  steps  the  smoothing  factor  actually  decreases,  so  the 
performance  of  the  solver  improves. 

2.5  Local  Mesh  Refinement 

The  resolution  is  increased  near  the  vortex  by  super¬ 
imposing  nested  uniform  grids  of  different  mesh  sizes,  in 
the  same  way  as  in  the  MUDBAR  model  (Fulton,  2000). 
Transfers  between  grids  are  more  complicated  due  to  the 
staggering  of  u  and  u.  The  only  conditions  applied  at  the 
grid  interfaces  are  the  same  as  the  open  boundary  condi¬ 
tions  applied  at  the  domain  boundaries,  with  the  specified 
quantities  (e.g.,  Vn  —  <p/c)  obtained  by  interpolating  in  time 
and  space  from  the  next  coarser  grid.  For  the  results  pre¬ 
sented  here  we  use  patches  of  fixed  size  which  are  moved 
to  follow  the  vortex  over  each  time  step. 


3.  NUMERICAL  RESULTS 

For  a  simple  test  case  we  use  a  small-scale  vortex 
(Gaussian  geopotential  with  e-folding  width  112  km)  em¬ 
bedded  the  sinusoidal  zonal  current  of  DeMaria  (1985). 
The  model  domain  is  a  square  of  side  length  4096  km  on  a 
/3-plane  at  20°  N  with  c  =  100  m/s,  using  a  base  grid  with 
/i  =  32  km  and  one  patch  with  h  =  16  km.  Figure  1  shows 
the  solution  (contours  of  geo  potential  and  wind  vectors)  at 
i  =  4  hours  and  i  =  48  hours;  the  smaller  square  shows  the 
boundary  of  the  grid  patch.  Since  the  initial  condition  is  in 
geostrophic  but  not  gradient  balance,  the  initial  unbalance 
generates  a  significant  gravity  wave  front  which  propagates 
cleanly  through  the  grid  interface  by  4  hours  and  out  of  the 
domain  long  before  48  hours.  There  is  no  evidence  of  false 
reflections  or  other  problems  at  the  grid  interface. 
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t=  4  hours 


Figure  1.  Sample  solution 
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